NUMERICAL STUDY OF THE SMALL DISPERSION LIMIT OF THE 
KORTEWEG-DE VRIES EQUATION AND ASYMPTOTIC SOLUTIONS 



T. GRAVA AND C. KLEIN 

Abstract. We study numerically the small dispersion limit for the Korteweg-de Vries (KdV) 
equation ut + Guux + ^Uxxx ~ for e ^ 1 and give a quantitative comparison of the numerical 
solution with various asymptotic formulae for small e in the whole (x, t)-plane. The matching of 
the asymptotic solutions is studied numerically. 



1. Introduction 

The behavior of solutions to Hamiltonian perturbations of hyperbolic and elliptic systems has 
seen a renewed interest in |28[ l3Uj . Specific integrable cases like the solution to the small 
dispersion limit of the Korteweg-de Vries (KdV) equation or the semiclassical limit of the nonlinear 
Schrodinger equation have been studied in detail in some part of the (x, t) plane in the seminal 
papers |501 [26l Wl\ . However some detail description of the solution in several critical regions of 
the (x, t) plane can be given for non-integrable Hamiltonian perturbations of hyperbolic or elliptic 
systems. It has been conjectured by Dubrovin and Dubrovin et al. 128\ [30l [29j that solutions can 
be approximated in one of the critical regimes by special solutions to the Painleve I equation and 
its hierarchy (see also [3]). In particular the universal nature of the critical behavior is remarkable. 
The conjecture has been rigorously proved in one specific case, that is for the Cauchy problem for 
the KdV equation with analytic initial data. Further critical behaviours have been observed in 
solutions to Hamiltonian perturbations of hyperbolic and elliptic equations (see e.g. [HdUE]). In 
particular in the Hamiltonian perturbations of hyperbolic systems, two other critical regimes have 
been observed: one of them is related to the second Painleve equation, while the other is solitonic 
because the local asymptotic behaviour is described by a train of solitons. The existence of such 
critical regimes has been rigorously proved for the Cauchy problem for the KdV equation. A review 
of these results as well as a new and improved numerical comparison of all asymptotic formulae 
with the numerical solution of KdV, is the subject of the present work. 

We consider the Cauchy problem for the KdV equation with small dispersion 

(1.1) ut + 6uux + e^uxxx=0, u{x,t = 0,e) = uo{x), e > 0, x G M, t G M+; 

uq{x) is real analytic negative initial data with sufficient decay at infinity and with a single negative 
hump (for detailed definition see [16|). For a much wider class of initial data than the one we 
consider, it is known that the KdV solution exists for all positive times t (see for example [23j). 
Up to a certain time tc the solution u{x, t, e) as e — )• can be approximated [50j by the solution to 
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the Cauchy problem for the (dispersionless) Hopf equation 

(1.2) Ut + Quux = 0, u{x,0) = uo{x), t > 0, 
which can be solved by using the method of characteristics in the form 

(1.3) u(x,t) = uo{0, x = 6tuo{0 + C- 
At time 

1 

tr. 



maxgg]R[-6no(^)] 

the Hopf equation reaches a point {xc, tc) of gradient catastrophe where the derivative of the Hopf 
solution blows up. Recently, it has been proved [Slj that for initial data uo{x) in the Sobolev space 
H'^ (M) , s > I , and for a larger class of equations than KdV 

(1.4) u{x,t,e) =u{x,t) + 0{e^), t < tc, 



where u{x,t) is the solution (1.3). The sub-leading term in the above expansion is determined 
explicitly. For t > tc, the Hopf solution (1.3) is multi- valued. However, the KdV solution is well- 
defined for all positive t and e > 0: the dispersive term e'^Uxxx regularizes the gradient catastrophe. 
For t slightly smaller than tc, the KdV solution u{x, t, e) starts to oscillate. For t > tc a zone of rapid 
modulated oscillations develops [131 [50]. In the (x,t)-plane, the oscillations take place in a cusp- 
shaped region (which depends on the initial data) as illustrated in Figjl] These oscillations in the 
limit e — )• are confined to a certain interval [x~ (t), x'^^t)], see Figure l2j The interval [x~ {t) , x^ {t)] 
is usually called Whitham zone because the oscillations are described inside this interval through 
the Whtiham equations ^0] (see below (2.7)). Furthermore the functions x^{t) and x^{t) are 
determined by the confiuent form of the Whitham equations (see subsections 2.3[ 2.4). 

Inside this cusp-shaped region, the exact one-phase solution to the KdV equation in terms of 
elliptic functions gives an asymptotic description of the oscillations, but on an elliptic surface where 
the branch points depend on x and t via the Whitham modulation equations as was proved by Lax 
and Levermore and Venakides 1 50\ [59] . This averaging procedure works well inside the Whitham 
zone, but has to be amended near the boundaries of the zone as was found numerically in (39] . 
Near the point of gradient catastrophe (xc, tc) , Dubrovin |28| conjectured that for a large class of 
equations containing KdV, the corresponding solution is asymptotically given in terms of a special 
solution of the second equation in the Painleve I hierarchy that we will call Pj equation. This was 
tested numerically in [31] and proven for KdV in [15] . Near the leading edge a multiscale expansion 
was presented for the oscillations in terms of a particular solution of the Painleve II equation in 
|40j . The validity of such an expansion has been proved rigorously in [161. Near the trailing edge 
[T7] gave an asymptotic solution in terms of a series of pulses of the shape of KdV solitons (see 
Fig I). 

It is remarkable that the KdV solutions can be approximated by Painleve equations in the critical 
regimes described above. Those Painleve equations have appeared in many branches of pure and 
applied mathematics during the past decades, for an overview we refer to j35j. 

The universal critical behaviour of KdV solutions should be seen in relation to the known univer- 
sality results in random matrix theory. For large unitary random matrix ensembles, local eigenvalue 
statistics turn out to be, to some extent, independent of the choice of the ensemble and indepen- 
dent of the reference point chosen |25[ 127] . Critical break-up times occur when the eigenvalues 
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Figure 1. Solution to the KdV equation for the initial data uq{x) = — sech^x and 
e = 10-^ 

move from a one-cut regime to a multi-cut regime. These transitions take place in the presence of 
singular points, of which three different types are distinguished [27]. Such singular points have the 
same singularity type which appears in the small dispersion limit of KdV. Singular interior points 
|34l [32| [GJ [T9| [20] show remarkable similarities with the leading edge of the oscillatory region for 
the KdV equation. The second possibility, which is related to the trailing edge for KdV, is that an 
interval in the spectrum shrinks and disappears afterwards. This case is also referred to as birth of 
a cut in unitary random matrix ensembles, ^ dU [SH [33]. The point of gradient catastrophe, i.e. 
the break-up point where the oscillations set in is comparable to a singular edge point in unitary 
random matrix ensembles. It was conjectured by Bowick and Brezin, and by Brezin, Marinari, and 
Parisi |9l [10] that local eigenvalue statistics in this regime should be given in terms of the Painleve I 
hierarchy. In |22j , it was proven that indeed double scaling limits of the local eigenvalue correlation 
kernel are given in terms of the Lax pair for the Pf equation. 

In this paper we implement numerically all known asymptotic approximations to the solution 
u{x, t, e) of the KdV equation as e — )• and test them quantitatively for a concrete example. With 
respect to previous works, the numerical methods are overall improved which allows to study a wider 
range of values for the small dispersion parameter e. The asymptotic formulae for the trailing edge 
are implemented for the first time as well as the terms of order e^/^ at the breakup point. We also 
address numerically the matching of different approximations. 

The paper is organized as follows: In sect. 2 we collect various asymptotic formulae for KdV 
solutions in the small dispersion limit. In sect. 3 we give a brief overview of the used numerical 
methods. In sect. 4 we study numerically the asymptotic description given by the one-phase KdV 
and the Hopf solution. At the leading edge, the multiscale solution in terms of a special solution 
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Figure 2. Whitham zone for the KdV equation for the initial data uo{x) = 
— sech^x. The shadow regions indicate the various asymptotic approximations in 
a neighbourhood of the Whitham zone. 

to the Painleve II equation is studied in sect. 5. In sect. 6 the same is done for the asymptotic 
solution at the trailing edge, and in sect. 7 for the point of gradient catastrophe. In sect. 8 we add 
some concluding remarks and an outlook on the connection formulae for the various asymptotic 
regimes. 

2. Asymptotic descriptions of the small dispersion limit 

In this section we summarize various asymptotic descriptions of the dispersive shocks which will 
be implemented numerically in the following. 

2.1. One-phase solution in the Whitham zone. Inside the cusp-shaped region in Fig. [2] for t 

slightly bigger than tc, the KdV solution u{x,t,e) can be approximated for small e, by the exact 
1-phase solution of KdV, where the branch points of the elliptic surface depends on x and t through 
Whitham equations. The one-phase solution of KdV can be written in terms of the Jacobi elliptic 
theta function and the complete elliptic integrals of the first and second kind E(s) and K(s) 
[501 EH US]: 

(2.1) ?x(x,t,e) ~/3i + /32 + /33 + 2a + 2e2— log??(J7(x,t);r), 
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where /3i > /32 > /^s, ^, a, and T have the form 



(2.2) 
(2.3) 



n{x,t) 



2eK{s) 
a{s) = -pi + (/3i 



/33) 



2t(/3i 
K{sy 



r = 



-f33)-q], 
K'js) 



/3i-/?3' 



Note that K'{s) = K{\/\ — s^), and is defined by the Fourier series 



,TTin T+2TTinz 



The formula for q in the phase in (2.2) is equal to j4m I26| 



(2.4) 



9(/3i,/32,/33) 



2^/2' 



vr 



dfidu 



where fiiu) is the inverse function of the decreasing part of the initial data uq. The formula (2.4) 



for q is valid as long as /33 does not reach the minimum value of the initial data uq. When 13^ 
reaches and goes beyond the negative hump it is also necessary to take into account the increasing 
part of the initial data /r, [39], [56] 



(2.5) 



1 

2^ 



dX 



+ 



1 



/32 V(/5i-A)(A-/32)(A-/33) 



Remark 2.1. In the formula (2.1 ) the term /3i + /32 + /33 + 2a is the weak limit of u{x, e) as e — )■ 



|50j and the term containing the 0- function describes the oscillations |26| I59j. The error term in 
the asymptotic expansion (2.1) should be of order 0(e). 



Formula (2.1) can be written also in terms of the Jacobi elliptic function dn in the form 

(2.6) u{x, t, e) ~ /32 + /33 - /3i + 2(/3i - l3s)dn\2K{s)Q + K{s)). 

For constant values of f3i > /32 > Ps, the right hand side of ( |2.6| ) is an exact solution of KdV. 
However in the description of the leading order asymptotics of u{x, t, e) as e — >• 0, the numbers 
Pi > P2 > Ps depend on x and t and evolve according to the Whitham equations [60] 



(2.7) 



da 9 ^ 
ot ox 



0, 



Pi + a 



+ 2{Pi + p2 + P3), i = 1,2,3, 



with a as in (2.3). 



The Whitham equations (2.7) can be integrated through the so-called hodograph transform, 



which generalizes the method of characteristics, and which gives the solution in the implicit form 

EI! 



(2i 



Vit + Wi 



where the Vi are defined in (2.7), and where w 



i = l,2,3, 
WiiPi, P2, Ps) for i 



1, 2, 3 is obtained from an 



algebro-geometric procedure by the formula [55] 
(2.9) w^ = Uv.,-2Y^Pk 

\ k=l J 



dq 

dpi 



1,2,3, 
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with q defined in (2.4) or (2.5). Solvability of (2.8) for initial data with a single hump was proved 
in 1561. 



Near the boundary of the oscillatory cusp-shaped region, neither the Hopf solution ( |l.3| nor the 
one-phase solution (2.1) gives a satisfactory description of the KdV solution u{x,t,e) as e — )■ 0. 
Three different transitional regimes can be distinguished: (1) the cusp point where the gradient 
catastrophe for the Hopf equation takes place and where f^i = (32 = = Uc, (2) the leading edge of 
the oscillatory zone where /32 = f^s, and (3) the trailing edge of the oscillatory zone where /3i = /32- 
We will illustrate in the next subsections that in the above three cases 



u{xc,tc) + 0{e7), 
u{x,t,e) = { u{x'{t),t) + 0{e^' 
uix+{t),t) + 0{l), 



near the point of gradient catastrophe (xc, tc) , 
near the leading edge, 
near the trailing edge, 



where u{x,t) is the solution of the Hopf equation and x^{t) are the boundaries of the Whitham 
zone. The sub-leading terms are described respectively by a Pj transcendent, a Pn transcendent, 
and a train of solitons. Clearly the above asymptotic expansions are not uniform in e. Connections 
formula need to be developed. 



2.2. Point of gradient catastrophe. It was conjectured in |28j and proved afterwards in |15j 
that the KdV solution u(x, t, e) as e — )■ near the point of gradient catastrophe (xc, tc) for the Hopf 
solution (1.3), is given in terms of a distinguished Painleve transcendent, namely a special smooth 
solution U {X, T) to the fourth order ODE 

(2.10) X = 6TU 



+ \ul + U Uxx + ^Uxxxx 



±oo. 

It is remarkable that 



This ODE is the second member of the Painleve I hierarchy, and we refer to it as the Pj equation. 
The relevant solution is real and has the asymptotic behavior 

(2.11) U{X,T) =^{\X\f/^ ^2T\X\-^/^ + 0{\X\-^), as X - 

for any fixed T G M, and has no poles for real values of X and T \21\ \28 
U{X, T) is an exact solution to the KdV equation 

(2.12) Ut + QUUx + Uxxx = Q- 

In a double scaling limit where e — t- and simultaneously x and t approach the point and the 
time of gradient catastrophe Xc and tc in such a way that the limits 

■ - Xc- 6Uc{t - tc) 



lim 
e 

6Urt — )• Xr 



6 

e7 



, lim 


'{t-tcY 

4 


, Uc — w(Xc, tc 


e 


. 67 





6Uctc t ^ tc 

exist and are bounded, the KdV solution has an expansion of the following form . Let 

- Xc - 6uc{t - tc) „ (t - tc) 



(2.13) 

with 
(2.14) 



X 



(A:)i/7e^ 



(fe)3/7e 



-f'/!inc)/6 
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and /^(n) the inverse of the decreasing part of the initial data. Then the solution of KdV is 
approximated by 



(2.15) 



u{x, t, e) 



Uc + 



2/7 



U{X,T) 



k) 63/"' (n,) 



1 



QUx + 2Uxx + 4?7" + 15T - 90T'Ux - 3XUUx - -XUxxx 



and Q{X,T) is the integral of U{X,T), 

Q = ^UxUxxx - Y^^lx + XU- 3TU^ + ^U^ + hjUl, Qx = U. 

2 

The correction term of order was rigorously derived using steepest descent analysis for the 
Riemann-Hilbert problem for KdV in |15| . Such an approximation was already discovered for the 
Gurevich-Pitaevskii solution of KdV in ISH. The correction of order e^/'' was derived in 1181. 



2.3. Leading edge. Near the leading edge at the left of the zone where the oscillations become 
small, a multiscale analysis and numerical results {40j showed that the envelope of the oscillations 
is asymptotically described by a particular solution to the second Painleve equation (P//) 

(2.16) q"{s) = sq + 2q^{s). 

The special solution we are interested in, is the Hastings-McLeod solution [44J which is uniquely 
determined by the boundary conditions 

(2.17) q{s) = x/-s/2(l +o(l)), ass^-oo, 

(2.18) g(s) = Ai(s)(l + o(l)), ass^+oo, 

where Ai(s) is the Airy function. Although any nonzero Painleve II solution has an infinite number 
of poles in the complex plane, it is known |44] that the Hastings-McLeod solution q{s) is smooth 
for all real values of s. 

The leading edge corresponds to the Whitham equations in the confluent case where 

Psix, t) = P2{x, t) = v{t), I3i{t) = u{t), 



see also (2.1). There exists a time t > tc such that for tc < t < t, the leading edge x (t) is 
determined uniquely by the system of equations [55l |l2] 

(2.19) x-it) = 6tu{t) + fLiuit)), 

(2.20) 6t + e{v{t);u{t)) = 0, 

(2.21) dj{vity,uit)) = 0, 
with u{t) > v{t) and with 



(2.22) 



9{v; u) = 0{v; u) 



1 



2V2J-1 



1 + m 1 

^ v + - 



m 



-u 



dm 



m 



Furthermore x~{t), u{t), and v{t) are smooth functions of t. Throughout the rest of the subsection, 
whenever we refer to u, we mean by this the solution of the system (2.19)-(2.21 ) for a given time 
tc < t < i, while we denote the solution of the KdV equation as u{x, t, e) and the solution of the 
Hopf equation as u{x, t). 
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The behaviour of the solution of KdV near the leading edge as e — )• is described as follows. 
Take a double scaling limit where e — )• and at the same time x — )• x~ (t) in such a way that 



(2.23) 



lim 

e ^ 

X — )• x~{t) 



X — X (t) 

I273 



remains bounded. In this double scaling limit, the solution u{x, t, e) of the KdV equation ( 1.1 ) with 
initial data uq has the asymptotic expansion |16j 



(2.24) u{x,t,e) 



u 



4eV3 
^1/3 



q (s(x, t, e)) cos 



G(x,t) 



+ 



X — X 



+ e^''^@i{x,t,e 
4e2/3 



6t + /[(n) c'^/^{u-v) 



{s{x,t,€)) sin 



Q{x,t] 



+ 0(6). 



Here x and v < u (each of them depending on t) solve the system (2.19), and the phase 0(x, t) is 
given by 



(2.25) 

Furthermore 
(2.26) 



@{x,t) = 2./^r^{x-x-) + 2 r{f^{^) + (St)^/i^di. 

J V 



52 

^/v^^^—-^9{v■,u)>{), s{x,t,e) 



X — X 



2/3' 



with 9 defined by (2.22), and q is the Hastings-McLeod solution to the Painleve II equation. The 
correction to the phase Qi{x,t,e) takes the form 



(2.27) Gi(x,t,e) 



.1/3 



+ p 



[v:u 



5p + 



6dl,9{v;u) 4{u-v) 

s{x,t,ef { dl,9{v;u) 



+ ■ 



+ 



2c^Ju — V 

3dl,9{v;u) 2{u-v) ' 6t + f'^{u) 
P = P{s) = -q^{s) - sq^{s) + q'{sf, 



where we used the notations 

Q = q{s), q' = q'{s), 

with s = s{x,t,e) and p'{s) = —q^- 

Remark 2.2. Note that the leading order term in the expansion ( |2.24 ) of u{x,t,e) is given by 
u{t) the solution of the Hopf equation at the leading edge. The second term in (2.24) is of order 
e^/^, while the remaining terms are of order e^/^. From the C(e^/^)-term, we observe that u{x, t, e) 
develops oscillations of wavelength 0{e) at the leading edge, the envelope of the oscillations is 

2 

proportional to the Hastings-McLeod solution q. If we let {x — x~{t))/e3 — )• —00 (so that x lies to 
the left of the leading edge) , the terms with the oscillations disappear due to the exponential decay 
of q, see (2.18). We are then left with only two terms in (2.24), which are the first two terms in 
the Taylor series of the Hopf solution u{x,t) near x~ . 
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Remark 2.3. The formula (2.24) can be obtained from a multiple scale analysis [40j from (2.1), 
letting 

2gi/3 2eV3 
Psix, t) = v{t) - —^q {s{x, t, e)) , ^2(2:, t) = v{t) + --^q {s{x, t, e)) , 



-.1/3 



/3i(t) = u(t) + 



X — X 



6t + 

2.4. Trailing edge. The trailing edge x~^{t) of the oscillatory interval is uniquely determined by 
the confluent form /3i = /32 = f and = u, v > u, of the equations (2.8), namely [S^ WI\ 

(2.28) x+{t) = 6tu{t) + fL{u{t)), 

(2.29) 6t + e{v{t);u{t)) = 0, 



(2.30) 
(2.31) 



v{t) 



u{t) 



(6i + 0{X; u{t)))^/X- u{t)dX = 0, 



where 9{v;u) has been defined in (2.22). 



Let x^ = x^{t), u = u{t), and v = v{t) solve the system (2.28)-(2.30), and let us take a double 
scaling limit where e — t- simultaneously with x — )• x^{t) in such a way that 



y 



2\A7" 



u ■ 



X — x^ 
elne 



remains bounded: there exists a real M > such that |y| < M. Then there exists i > tc such that 
for tc < t < t, we have the following expansion for the KdV solution u{x,t,€) in the double scaling 
limit il7j, 

[Ml 

(2.32) u{x,t,e) = u + 2{v-u)^sech'^{Xj) + 0{eln'^e), 

j=0 

where \M~\ is the smallest integer > M, 



(2.33) 



hi 



^(^ - 2/ + j) Ine - HVl^hj) + log 7, 



2i 



7r4VJ 



5 . 

, 'J = 4:{v - u)4y'-djj0{v;u), 



and 9 is given by (2.22). Observe that hj are the normalization constants of the Hermite polyno- 
mials. 



Remark 2.4. Observe that each term in the sum of ( 2.32 ) generates a pulse with amplitude 2{v—u) 
for y near a half positive integer which can be seen as a soliton. Indeed the term sech^ [Xj ) is of the 

order 0(1) for y = j + 1/2. For y = j or j + 1, it already decreased to order 0(e2). For y = j — \ 
or j + |, the contribution of sech^(Xj) is absorbed by the error term 0(eln^ e). Clearly, since j is 
nonnegative, the solitons appear only for y positive, that is for a; < x"*", namely inside the Whitham 
zone. 
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Remark 2.5. The phase defined m (2.2) satisfies the formal limit 

and from the Whitham solution (2.8) one obtains in the limit /3i — )■ /32 

'/3i-/32 



X — X (t) ~ 



--dj{v;u) ( 



log 



The Jacobi elliptic function dn — )• sech as the modulus s — )• 1 and 



K{s) ~ ^ log 



1 



as s 



1. 



Therefore the formal limit of the solution (2.6) in terms of elliptic functions as s 
u, gives 



1, ^1,^2 ->■ V, 



(2.34) 



u{x, t, e) ~ ti + 2{v — u) sech^ 

he perio 

'/3i-/32 



x — X' 



1, 



-Vv-u + {j + -)log 



1 



for any positive integer j, due to the periodicity of the elliptic function dn. Choosing 

46 



dy9{v, u)y/v — u 



and inserting it in (2.34) one can partially reproduce the formula (2.32) in the sense that all 
the terms in the phase Xj defined in (2.33) can be reproduced except the one containing the 
normalization constants of the Hemite polynomials hj. We would like to remark that the formal 



limit (2.34) has appeared several times in the literature, but such limit does not describe the small 



dispersion solution of KdV near the trailing edge, since the limiting value of the phase (2.2) does 
not give the right result. 

3. Numerical Methods 

The numerical task in treating the small dispersion limit of KdV and various asymptotic formulas 
consists in solving the KdV equation itself, certain ODEs of Painleve type for a given asymptotic 
behavior, and of the Whitham equations for which the implicit solution (2.7) exists. We will 



summarize in this section how these different tasks are solved numerically, and how we control the 
numerical accuracy. 

3.1. KdV solution. Since critical phenomena are generally believed to be independent of the 
chosen boundary conditions, we study a periodic setting in the following. This also includes rapidly 
decreasing functions which can be periodically continued as smooth functions within the finite 
numerical precision. This allows to approximate the spatial dependence via truncated Fourier series 
which leads for the studied equations to large stiff systems of ordinary differential equations (ODEs), 
see below. The use of Fourier methods not only gives spectral accuracy in the spatial coordinates 
(the numerical error in approximating smooth functions decreases faster than any power of the 
number A'^ of Fourier modes) , but also minimizes the introduction of numerical dissipation which is 
important in the study of the purely dispersive effects we are interested in here. In Fourier space. 



equation (1.1) has the form 
(3.1) 



vt 



Lv + 'N{v,t), 
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where v denotes the (discrete) Fourier transform of u, and where L and N denote hnear and non- 
hnear operators, respectively. The resulting system of ODEs consists in this case of stiff equations 
where the stiffness is related to the linear part L (it is a consequence of the distribution of the 
eigenvalues of L), whereas the nonlinear part contains only low order derivatives. In the small 
dispersion limit, this stiffness is still present despite the small term in L. This is due to the fact 
that the smaller e is, the higher wavenumbers are needed to resolve the rapid oscillations. 

Loosely speaking a stiff system is a system for which explicit numerical schemes as explicit Runge- 
Kutta methods are inefficient, since prohibitively small time steps have to be chosen to control 
exponentially growing terms. The standard remedy for this is to use stable implicit schemes, which 
require, however, the iterative solution of a system of nonlinear equations at each time step which 
is computationally expensive. In addition the iteration often introduces numerical errors in the 
Fourier coefficients. Thus we used in |39j an integrating factor method, where the linear stiff part 
is explicitly integrated. This can be conveniently done here since the operator L corresponding 
to the third derivative with respect to x is diagonal in Fourier space. As was shown in j45] . 
integrating factor methods can suffer from order reductions, which means that the actual decrease 
of the numerical error with the numerical resolution is much lower than the classical order of the 
used method. This was confirmed for the small dispersion limit of KdV in |48j . There it was also 
shown that exponential time differencing (ETD) schemes are very efficient for KdV. ETD schemes 
were developed originally by Certaine [12j in the 60s, see [46j for a comprehensive review. The 



basic idea is to use equidistant time steps h and to integrate equation (3.1) exactly between the 



time steps i„ and tn+i with respect to t. With v{tn) = Vn and v(t„+i) = fn+i, we get 



Vn+l = e^^Vn 



+ / e^'^^-^^^{v{tn + T),tn + T)dT. 

Jo 



The integral will be computed in an approximate way for which different schemes exist. We use 
here a Runge-Kutta method of classical order 4 due to Cox-Matthews [2l]. As in [T| for the 
Camassa-Holm equation, this approach could be amended by identifying two regimes t G [0, ti] and 
t G [ti,tend] with ti <^ tc- A much larger time step can be used in the first regime than in the 
second where the rapid modulated oscillations appear. We do not use this approach here since it 
was not necessary for the considered values of e. 

The accuracy of the numerical solution is controlled via the numerically computed conserved 
energy of the solution 

(3.2) E[u] = / {2u^ - e^ul)dx, 

which is an exactly conserved quantity for KdV. Numerically the energy E will be a function of 
time. We define AE := \{E{t) — £'(0))/£'(0)|. It was shown in ^48j that this quantity can be used 
as an indicator of the numerical accuracy if sufficient resolution in space is provided. The quantity 
AE typically overestimates the precision by two orders of magnitude. Since the numerical error 
has to be clearly smaller than the difference between KdV solution and the asymptotic descriptions 
we want to test (which give at best descriptions of order e) we are interested in a numerical value 
clearly below the smallest considered value of e. To ensure this we will always ensure that the 
modulus of the Fourier coefficients of the final state decreases well below 10~^ (thus providing the 
needed resolution), and that the quantity AE is smaller than 10~^ (in general it is of the order of 
machine precision, i.e. 10~^^). 
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We consider in the following always the example uo = — sech^x and values of 

e = 10-\ 10-1-25,... 10-3-^ 

For the smallest values of e, we use = 2^^ Fourier modes and A^^ = 4 * 10^ time steps; for larger 
values of e between 2^^ and 2^'' Fourier modes, and between 10'^ to 10^ time steps. The oscillatory 
zone for e = IQ-^'^ can be seen in Fig. [sj The Fourier coefficients for this solution are shown in 
Fig-i 




-3.2 -3 -2.8 -2.6 -2.4 -2.2 -2 



Figure 3. The oscillatory zone in the solution to the KdV equation for the initial 
data uo{x) = — sech^x and e = IQ-^'^ for t = 0.4. The oscillations are so rapid, that 
they are graphically difficult to represent though they are numerically well resolved. 



3.2. Numerical solution of the Whitham equations and of the Hopf equation. The 

Whitham equations (2.7) are solved for given initial data by inverting the hodograph transform 
(2.8) to obtain (3i > (32 > as a function of x and t, and similarly for the implicit solution of the 
Hopf equation (1.3). Since the hodograph transform becomes degenerate at the leading and trail- 
ing edge we solve the system (2.19)-(2.21 ) and (2.28)-(2.30) instead of (2.8) to avoid convergence 
problems. 

These equations are of the form 

(3.3) Si{{yi},x,t) = 0, i = l,...,M, 

where the Si denote some given real function of the yi and x, t. The task is to determine the yi in 
dependence of x and t. To this end we determine the yt for given x and t as the zeros of the function 
S := X^i^i Sf. This is done numerically by using the algorithm of [39] which is implemented as the 
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2 3 

X 10' 



Figure 4. Fourier coefficients for tlie solution in Fig. [sj 



function fminsearch in Matlab. Tlie algoritlim provides an iterative approacli wliicli converges in 
our case rapidly if the starting values are close enough to the solution (see below how the starting 
values are chosen). We calculate the zeros to the order of machine precision. 

For a given t > tc, we always first solve the system (2.19)-(2.21 ) to obtain the leading edge 
coordinate x~{t) and 

Similarly we solve the equations (2.32 ) for and = and which fixes the interval [x~ , x~^] . 
This interval is subdivided into a number of points x„, n = 1, . . . , N^- In contrast to [39], we choose 
the Xn to be related to Chebyshev collocation points Ij = cos{jTT/Nc), j = 0,1 . . . , Nc, to allow for 
better interpolation formulas. Since the polynomial interpolation we will use works best for smooth 
functions, we use the analytic knowledge that /32 ~ /?3 ~ \/ x — x^{t) for x ~ and similarly 

/3i ~ /32 ~ Y^x+(t) — X for X ~ x~^{t). Thus we put for j = 0, . . . , Nc 



x 



and 



-it) + 



~{t)-x-{t)il + l,f 



-it)-x~it){i-i,r^ 



1 



x(.[x-^t),-{x-it)+x+m 



^ e ilix'it) +x+{t)),x+{t)]. 



2 4 

For given Xj and t, the Whitham equations are solved as discussed in [HD]. Thus the /3j are sampled 
on Chebyshev collocation points which can be used to obtain an expansion of these functions in 
terms of Chebyshev polynomials, see for instance [SB]- As for Fourier series, the order of magnitude 
of the modulus of the coefficient of the highest order polynomial gives for smooth functions an 
indication of the numerical resolution. For our example the Chebyshev coefficients decrease well 
below 10~^ with Nc = 64 which is more than sufficient for our purposes. To obtain machine 



precision, the integrals in (2.8) would have to be computed as described in |39] with higher precision 
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Figure 5. Solutions of the Whitham equations (2.8) for the initial data uq = 
— sech^x for t = 0.4. 

At intermediate values of x G [x~{t),x~^{t)], the /3j are obtained from the values on the collocation 
points via numerically stable barycentric Lagrange interpolation, see |2j, which is essentially an 
efficient implementation of the Lagrange polynomial for Chebyshev collocation points. 

3.3. Painleve transcendents. The asymptotic solutions near the breakup point and the leading 
edge are given by pole-free solutions with a given asymptotic behaviour for x — t- itcxD to the 
equation and the Painleve II equation respectively. A way to solve these equations is to give a series 
solution to the respective equation with the imposed asymptotics that is generally divergent. These 
divergent series are truncated at finite values x, xi < Xr at the first term that is of the order of 
machine precision. The sum of this truncated series at these points is then used as boundary data, 
and similarly for derivatives at these points. Thus the problem is translated to a boundary value 
problem on the finite interval [xi , Xr] ■ 

In [lO] we used for the P| solution a collocation method with cubic splines distributed as bvp4 
with Matlab, in [H] for the Hastings-McLeod solution of P// a Chebyshev collocation method with 
a fixed point iteration. Here we use again a Chebyshev collocation method for both equations. As 
for the Whitham equations above, the solution of the ODEs is sampled on Chebyshev collocation 
points Xj, j = 0, . . . , Nc which can be related to an expansion of the solution in terms of Chebyshev 
polynomials. Since the derivative of a Chebyshev polynomial can be again expressed in terms of a 
linear combination of Chebyshev polynomials, the action of the derivative operator on the Hilbert 
space of Chebyshev polynomials is equivalent to the action of a matrix on this space. This leads 
to the well known Chebyshev differentiation matrices, see for instance j57]. Thus for the numerical 
solution, in an ODE of the form F{u, dxU, ..) = 0, n is replaced by the vector u{xj), j = 0, . . . , Nc 
and dx by the differentiation matrix. The ODE is in this setting replaced by A^c + 1 algebraic 
equations. The boundary data are included via a so-called r-method: The equations for j = 
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and for j = (for the fourth order equation j = 0, 1, A^c — li-^c) are replaced by the boundary 
conditions. The resulting system of algebraic equations is solved with a standard Newton method. 
The convergence of the solutions is in general very fast. We always stop the Newton iteration when 
machine precision is reached. Again the highest Chebyshev coefficients are taken as an indication of 
sufficient resolution of the solutions (they have to reach machine precision) . A similar approach had 
been used in [7] for the Hastings-McLeod solution. Solutions to the Painleve II equation have been 
computed as the solution of a Riemann-Hilbert problem in |53j . Certain Painleve transcendents 
can be expressed in terms of Predholm determinants which can be computed with the methods of 
[8]. For the study of Painleve solutions with poles in the complex plane, an approach based on 
Pade approximants has been presented in [37] . 

The Hastings-McLeod solution and the special solution to the Pj equation for various values of 
t can be seen in Fig. [6} 




Figure 6. Hastings-McLeod solution of the Painleve II equation on the left, and 
the special solution to the Painleve 12 equation for several values of t on the right. 



4. Numerical solution of KdV, Whitham and Hope equations 



In this and the following sections, we will numerically solve the KdV equation (1.1) for the 
initial data uq{x) = — sech^x for the values e = 10^^, 10~^'^^, . . . , 10~'^'^, and compute the various 
asymptotic descriptions of the small dispersion limit from sect. 2. We will study the validity of 
these asymptotic descriptions in various regions of the (x, t)-plane. To obtain the e-dependence of 
a certain quantity A, we perform a linear regression analysis for the dependence of the logarithms, 
In ^ = a In e + 6. This allows to obtain numerically the scaling of the difference between numerical 
and asymptotic solutions also for the cases where no analytic behavior is yet known. We first 
consider the asymptotic description based on the Hopf solution outside the Whitham zone and the 
one-phase KdV solution inside the zone with branch po of the elliptic surface given by the Whitham 
equations. 

Outside the Whitham zone, the Hopf solution for the same initial data as the KdV solution gives 
an asymptotic description of the latter. Inside the Whitham zone, the one-phase KdV solution 
provides an asymptotic solution. 



16 



T. GRAVA AND C. KLEIN 



Before breakup: For times much smaller than the critical time, we find that the L^o norm of the 
difference between Hopf and KdV solutions decreases as e"^. More precisely we find by linear 
regression an exponent a = 1.9987 with correlation coefficient r = 0.999995 and standard deviation 
a = 0.0051. 

At breakup, t = tc-' For times close to the breakup time, the Hopf solution develops a gradient 
catastrophe. The largest difference between Hopf and KdV solution can be found close to the 
breakup point. We determine the scaling of the Loo norm of the difference between Hopf and KdV 
solutions on the whole interval of computation. We find that its scaling is compatible with e^/^ 
as conjectured in [28] and proven in [15]. More precisely we find in a linear regression analysis 
a = 0.2929 (2/7 = 0.2857. . .) with a correlation coefficient r = 0.99996 and standard deviation 
a = 0.0022. 

After breakup: For times much greater than the critical time of the Hopf solution, we find that 
the asymptotic solution given by Hopf solution and the one-phase KdV solution via the Whitham 
equations gives a very good description of the KdV solution. Thus it is necessary to plot the 
difference between these solutions as done in Fig. [7] 

It can be seen that this difference is not uniform in x, and that this applies also for the decrease 
with e. The approximation is very good close to the centre of the Whitham zone, but much worse at 
the edges. We define the interior part of the zone tentatively as the interval symmetric to the centre 
of half the length of the zone. The results do not depend on whether this zone is taken slightly 
smaller or bigger. We find that the Lqo norm of the difference of the numerical KdV solution 



and the one-phase KdV solution (2.1) decreases there roughly as e. More precisely we find in a 
linear regression analysis a = 0.98 with a correlation coefficient r = 0.998 and standard deviation 
a = 0.047. 

At the leading edge we find that the error is always biggest close to the boundary of the Whitham 
zone. In an interval symmetric to this boundary with the same length as the above interior zone, 
we find that the diff eren ce between KdV solution and asymptotic solutions via Hopf and the one- 



phase KdV solution (2.1 ) scale roughly as e^/^. More precisely we find in a linear regression analysis 
a = 0.33 with a correlation coefficient r = 0.999 and standard deviation a = 0.012. 

The situation at the trailing edge is more complicated. It can be seen in Fig. [7] that the difference 



between the KdV solution and one-phase KdV solution (2.1 ) is almost constant (roughly 0.04) there. 
Notice that this error of order 0{1) which is supposed to appear in a zone of width elne close to 
the trailing edge of the Whitham zone was not seen in [39] because of a lack of resolution. This 
is one of the reasons why we redid the computations with a considerably higher resolution. It can 
also be seen in Fig. [7]that the 0(1) oscillation is moving closer and closer to the edge with smaller 
e as expected. The difference between KdV and Hopf solutions close to the trailing edge decreases, 
however, roughly as ^/e. More precisely we find in a linear regression analysis a = 0.54 with a 
correlation coefficient r = 0.997 and standard deviation a = 0.03. 

5. Leading edge 



In this section we study numerically the asymptotic formula (2.24) via the Hastings-McLeod 
solution which approximates the KdV solution at the leading edge as e — )• 0. We will refer to this 
asymptotic solution as P// asymptotics. We identify the zone, where the P// asymptotics gives a 



better description of KdV than the Hopf (1.4) or the one-phase KdV solution (2.1) and study the 
e-dependence of the errors. 
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Figure 7. The blue line describes the difference between the numerical solution 



of the KdV equation and the asymptotic formula (2.1) for the initial data uq{x) = 
— 1/ cosh^x and for t = 0.4. The green lines represent the difference between the 



numerical solution of the KdV equation and the Hopf solution (1.3) 



In Fig. [S] we show the KdV solution, the asymptotic solution via Whitham and Hopf and the 
P//- asymptotics near the leading edge of the Whitham zone. It can be seen that the one-phase 
KdV solution gives a very good description in the interior of the Whitham zone as discussed above, 
whereas the P// asymptotics gives as expected a better description near the leading edge. 

In Fig. [ojthe KdV solution and the P// asymptotics are shown in one plot for e = 10~^. It can 
be seen that the agreement near the edge of the Whitham zone is so good that one has to study 
the difference of the solutions. The solution only gives locally an asymptotic description and is 
quickly out of phase for larger distances from the leading edge, whereas the amplitude is roughly 
of the right size. 
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Figure 8. The figure shows in the upper part the numerical solution to the KdV 
equation for the initial datum uq = — sech^x and e = 10~^ at t = 0.4, in the middle 
the corresponding asymptotic solution in terms of Hopf and one-phase KdV solution, 



and in the lower part the P// asymptotic solution (2.24). 
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Figure 9. The numerical solution to the KdV equation for the initial datum uq = 
— sech^x and e = 10~^ at t = 0.4 in blue and the corresponding P// asymptotic 



solution (2.24) in green. 
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Figure 10. The difference between the numerical solution to the KdV equation for 
the in itial datum uo = — sech^x at t = 0.4 and the corresponding multiscale solution 
(2.24) for four values of e. Notice the scaling of the x and A axes with a factor e^/^ 



and e respectively to take care of the expected scalings in x and of the shown error 
next to the leading edge of the Whitham zone. 



The difference between KdV solution and the P// asymptotics is shown for several values of e in 



Fig. 10 It can be seen that the error close to the Whitham edge is almost constant. The scales 
in X and A are rescaled by a factor e^/'^ and e respectively which is the expected scaling behavior 
of the zone, where the multiscale solution should be applicable, and of the expected error. It can 
be seen that with these rescalings the error is of the same order for different values of e. A linear 
regression analysis for the logarithm of the difference A between KdV and multiscale solution in 
the interval [x~ — e^/'^, x~ + e^/'^] gives a scaling of the form A oc e"^ with a = 1.00 with standard 
deviation aa = 0.004 and correlation coefficient r = 0.99999. The result is almost the same in a 
larger interval, e.g., [x~ — 2e'^^^,x~ +2e^/'^] with just a slightly worse correlation. The found scaling 
is thus as expected of order e. 
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Figure 11. In Fig. (a) the difference between the numerical solution to the KdV 
equation for the initial datum uq = — sech^x at t = 0.4 for e = 10"^ and the 



corresponding P// asymptotic solution (2.24) in blue, and the difference between 



KdV and Hopf and one-phase KdV solution in green. In Fig. (b) the edges of 
the zone where the P// asymptotic solution (2.24) provides a better asymptotic 



description of KdV than the Hopf or the one-phase KdV solution in dependence of 
e. 



As can be already seen from Fig. |8} the multiscale solution gives a better asymptotic description 
of KdV near the leading edge of the Whitham zone than the Hopf and the one-phase KdV solution. 



This is even more obvious in Fig. 11a where the difference between KdV and the asymptotic 
solutions is shown. 

This suggests to identify the regions where each of the asymptotic solutions gives a better de- 
scription of KdV than the other. The results of this analysis can be seen in Fig. 12 This matching 
procedure clearly improves the KdV description near the leading edge. We also show the difference 
between this matched asymptotic solution and the KdV solution for two values of e. Visibly the 
zone, where the solutions are matched, decreases with e. 

There is a certain ambiguity in the precise definition of this matching zone due to the oscillatory 
character of the solutions. The limits of the matching zone for several values of e can be seen in 



Fig. lib Due to the lower number of oscillations in the Hopf region, the matching zone extends 



much further into this region than in the Whitham region. There does not appear to be a clear 
scaling law for the width of this zone. It can be already seen in Fig. 12 that the error at the 
matching does not scale with e. In fact we find a scaling close to with o ~ 2/3 (in the Whitham 
zone we find a = 0.63 and aa = 0.015 with r = 0.9995, and in the P// zone a = 0.60 and aa = 0.063 
with r = 0.99). Thus it is not possible to obtain an error of order e up to the trailing edge. It is 
clear that analytic connection formulae between the two asymptotic solutions must be established 
to obtain an error of order e in the shown range. 
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Figure 12. In the upper part of Fig. (a) one can see the difference between the 
numerical solution to the KdV equation for the initial datum uq = — sech^x and 
e = 10~^ at t = 0.4 and the corresponding asymptotic solution in terms of Hopf 
and one-phase KdV solutions. The lower part shows the same difference, which is 
replaced close to the leading edge of the Whitham zone by the difference between 
KdV solution and the P// asymptotic solution ( 2.24[ ) (shown in red where the error 
is smaller than the one shown above) . The figures in (b) show the same situation as 
in the lower part of (a) for two values of e. Notice the rescaling of the A axis with 
a factor e, the expected scaling of the error. 



6. Trailing edge 



In this section we study numerically for times greater than the critical time the soliton asymptotic 
formula (2.32) that approximates as e — >• the solution of KdV near the trailing edge of the 



oscillatory zone. We identify the zone, where this asymptotic formula gives a better description of 



KdV than the one-phase KdV (2.1) and Hopf (1.4) solutions and study the e-dependence of the 
errors. 



In Fig. 13 we show the KdV solution, the asymptotic solution via Whitham and Hopf and the 
soliton asymptotics near the trailing edge of the Whitham zone. As before the one-phase KdV 
solution gives a very good description in the interior of the Whitham zone, whereas the soliton 
asymptotic formula gives as expected a better description near the trailing edge. 
In Fig. 
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the KdV and the multiscale solution are shown in one 10"^ It 

seen that the agreement very close to the boundary of the Whitham zone is once more so good 
that the difference of the solutions has to be studied. The solution only gives locally an asymptotic 
description, and the quality of the approximation is not symmetric around the critical point. 
The difference between KdV solution and the soliton asymptotic solution is shown for several 



values of e in Fig. 15 The scales in x and A are both rescaled by a factor e which is the expected 
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Figure 13. The figure shows in the upper part the numerical solution to the KdV 
equation for the initial datum uq = — sech^x and e = 10~^ at t = 0.4, in the middle 
the corresponding asymptotic solution in terms of Hopf and one-phase KdV solution, 



and in the lower part the multiscale solution (2.32). 



scaling behavior of the zone (numerically e and elne are indistinguishable), where the soliton 
asymptotic solution should be applicable, and of the expected error. It can be seen that with these 
rescalings the error is of the same order for different values of e. A linear regression analysis for the 
logarithm of the difference A between KdV and multiscale solution in the interval [x"*" + e In e, x"*" — 
elne] gives a scaling of the form A oc e" with a = 1.07 with standard deviation 0"^ = 0.056 and 
correlation coefficient r = 0.998. The found scaling is thus compatible with elne. 



It can be seen in Fig. 16a where the difference between KdV and the asymptotic solutions is 
shown, that soliton asymptotic solution gives a much a better description of KdV near the trailing 
edge of the Whitham zone than the Hopf and the one-phase KdV solution. 

Again we can identify the regions where each of the asymptotic solutions gives a better description 
of KdV than the other. The results of this analysis can be seen in Fig. 17a This matching procedure 



clearly improves the KdV description near the trailing edge. In Fig. 17b we see the difference 
between this matched asymptotic solution and the KdV solution for two values of e. Visibly the 
zone, where the solutions are matched, decreases with e. 

Once more there is no precise definition of this matching zone due to the oscillatory character of 
the solutions. We determine it as the point where the curves of the differences intersect, or where 
they come closest, before one error dominates the other for all smaller respectively larger values 
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of X. The limits of the matching zone for several values of e can be seen in Fig. |16b There does 



not appear to be a clear scaling law for the width of this zone. It can be already seen in Fig. 17b 



that the error in the matching zone does not scale with e as close to the boundary of the Whitham 
zone. In fact we find a scaling close to e^^^ in both cases, but the correlation is not very good. As 
for the leading edge, it is necessary to establish analytical connection formulae. 



7. Point of gradient catastrophe 



as e 



In this section we study numerically the approximation (2.15) to the solution u{x,t,e) of KdV 
near the point of gradient catastrophe (xc, tc) for the solution of the Hopf equation. We 



identify the zone, where the P| asymptotic formula (2.15) gives a better asymptotic description of 
KdV than the Hopf or the one-phase KdV solution and study the e-dependence of the errors. We 
qualitatively study for a time t > tc close to tc how the various multiscale approximations perform. 



7.1. Critical time. For the initial datum uo{x) = — sech^x the critical time is t 
and the critical point Xc = — \/3/2 + ln((-v/3 — l)/\/2) ■ 
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. = V3/8 ~ 0.2165 
we show the KdV 



-1.5245. In Fig. 

solution, the Hopf solution and the multiscale solution near the critical point of the Hopf solution 
at the critical time. As before the Hopf solution gives a very good description for \x — Xc\ ^» 0, 
whereas the Pj asymptotic solution gives as expected a better description near the critical point. 
The following figures are always symmetric with respect to x, 
In Fig. 
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the KdV solution and the P 



J asymptotic solution are shown in one plot for e 



10 



-2 



It can be seen that the agreement very close to the critical point of the Hopf solution is again 
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Figure 1 5 . The difference between the numerical solution to the KdV equation for 
the initial datum uo = —sech^x at t = 0.4 and the corresponding soliton asymptotic 



solution (2.32) for four values of e. Note the scaling of the x and A axes with a 



factor e to take care of the expected scaling in x and of the shown error next to the 
trailing edge of the Whitham zone. 



so good that the difference of the solutions has to be studied. The solution only gives locally an 
asymptotic description. 

The difference between KdV solution and Pj asymptotic solution is shown for several values of 
e in Fig. 20 The scales in x are rescaled by a factor e^/'', and the ones for A by a factor e^/^ 



respectively which is the expected scaling behavior of the zone, where the Pj asymptotic solution 
should be applicable, and of the expected error. It can be seen that with these rescalings the error 
is of the same order for different values of e, at le ast cl ose to the critical point. We also show in 
this figure the different behaviour of the terms in (2.15) in order e^/'' and order e"^/^. The former 
is not symmetric with respect to the critical point. In fact the approximation is better on the side 
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Figure 16. In Fig. (a), the difference between tlie numerical solution to the KdV 
equation for the initial datum uq 



— sech X and e = 10~ at t = 0.4 and the 
corresponding soliton asymptotic (2.32) in green, and the difference between KdV 
and Hopf and one-phase KdV solution in blue. In Fig. (b) the edges of the zone 



where the soliton asymptotic (2.32 ) provides a better asymptotic description of KdV 



than the Hopf or the one-phase KdV solution in dependence of e. 



where the oscillations appear. However if one studies positive initial data as in [31], the oscillations 
are on the side of the critical point where the approximation in terms of the P| solution is worse. 

As is clear from Fig. 18, the multiscale solution gives a better asymptotic description of KdV 
near the critical point than the Hopf. This is even more obvious in Fig. 21a| where the difference 
between KdV and the Hopf solution is shown. 

Again we can identify the regions where each of the asymptotic solutions gives a better description 
of KdV than the other. The results of this analysis can be seen in Fig. 22a This matching procedure 



clearly improves the KdV description near the critical point. In Fig. 22b we see the difference 
between this matched asymptotic solution and the KdV solution for two values of e. Visibly the 
zone, where the solutions are matched, decreases with e (notice the rescaling of the axes with e). 

Once more there is no precise definition of this matching zone due to the oscillatory character of 
the solutions. We choose it as in the previous sections as given where the curves of the differences 
intersect, or where they come closest before one error dominates the other. The limits of the 
matching zone for several values of e can be seen in Fig. 21b There does not appear to be a 



clear scaling law for the width of this zone. A linear regression analysis for the logarithm of the 
difference A between KdV and multiscale solution in the matching zone gives a scaling of the form 
A oc with a = 0.586 (4/7 ~ 0.5714) with standard deviation aa = 0.06 and correlation coefficient 
r = 0.99 for the terms up to order e^/^ and with a = 0.62 (5/7 ~ 0.7143) with standard deviation 
o"a = 0.09 and correlation coefficient r = 0.98 for the terms up to order e^/^. The found scaling is 
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(A) (B) 

Figure 17. In the upper part of Fig. (a) one can see the difference between the 
numerical solution to the KdV equation for the initial datum uq = — sech^x and 
e = 10~^ at t = 0.4 and the corresponding asymptotic solution in terms of Hopf 
and one-phase KdV solution. The lower part shows the same difference, which is 
replaced close to the trailing edge of the Whitham zone by the difference between 
KdV solution and the soliton asymptotic ( |2.32| ) (shown in red where the error is 
smaller than the one shown above). In Fig. (b) the same situation as in the lower 
part of (a) is shown for two values of e: 10"^, 10~^. The A-axis is rescaled by a 
factor e. 



thus compatible with the expected e^/^ and e^/^ respectively. The error in the Hopf zone at the 
limits of the matching zone are of the same order. As before it would be interesting to study the 
connection formulae between the Hopf and P| zone. 

7.2. Close to the critical time. It is an interesting question to study how the various multiscale 
approximations perform for a time t greater than the critical time with t ^ tc, where all previously 
studied asymptotic formulae should be applicable. We consider just one value of e (e = 10~^) since 
the various multiscale expansions use different e-dependent rescalings of the time. We consider the 
time t = 0.23 > tc ~ 0.2165. 

First we study the situation in the vicinity of the leading edge (x = —1.6051 . . .). In Fig. 23 



the KdV solution can be seen in this case as well as the asymptotic solution in terms of Hopf and 
one-phase KdV solution, the multiscale solution (2.24) near the leading edge and the multiscale 



solution (2.15) close to the critical point. 



In Fig. 24 the same solutions can be seen in one figure. 

It can be seen already from these figures or from the plot of the differences between the asymptotic 



solutions and the KdV solution in Fig. 24b that the asymptotic solution in terms of Hopf and one- 
phase KdV solution performs worst close to the leading edge, and that the multiscale solution 
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Figure 18. The figure shows in the upper part the numerical solution to the KdV 
equation for the initial datum no = — sech^x and e = 10~^ at t = tc, in the middle 
the corresponding Hopf solution, and in the lower part the P| asymptotics (2.15). 



(2.15) in terms of the P| transcendent is most satisfactory. The multiscale solution (2.24) is only 
very close to the leading better than the asymptotics. It captures qualitatively the oscillations 
in the Hopf region, but is not oscillating around the Hopf solution as KdV. As can be seen from 
(2.24) the reason for this is that the amplitude of this solution is divided hy u — v which tends to 
at the critical point. The P| asymptotics also quickly becomes out of phase in the Hopf region. 
Thus to obtain a satisfactory description of the small dispersion limit of KdV close to the critical 
point, one has to study connection formulae between the various asymptotic solutions. 

The situation near the trailing edge {x = —1.5757. . .) is similar. In Fig. 25 the KdV solution 



can be seen in this case as well as the asymptotic solution in terms of Hopf and one-phase KdV 
solution, the multiscale solution (2.32) near the leading edge and the multiscale solution (2.15) 
close to the critical point. 

In Fig. 26 the same solutions can be seen in one figure. These figures as well as the plot of the 
differences between the asymptotic solutions and the KdV solution in Fig. 26b show that the as- 
ymptotic solution in terms of Hopf and one-phase KdV solution performs worst close to the trailing 
edge, and that the P| asymptotic solution (2.15) performs best. The soliton asymptotic solution 
(2.24) is better than the P| asymptotics only very close to the trailing edge. The P| asymptotics 
also quickly becomes unsatisfactory in the Hopf region. Thus to obtain a better description of 
the small dispersion limit of KdV close to the critical point, one has to study connection formulae 
between the various asymptotic solutions. 
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Figure 19. The numerical solution to the KdV equation for the initial datum 
uq = — sech^x and e = 10~^ at t = tc in blue and the corresponding asymptotic 
solution in green. 



8. Outlook 



The numerical results of the previous sections have shown that the proposed asymptotic descrip- 
tions lead in fact to an error of order e in various regions of the x, t-plane where the respective 
formulae are supposed to hold. However the asymptotic description at the trailing edge of the os- 
cillatory zone and the point of gradient catastrophe are characterized by errors of higher order. In 
order to obtain a complete analytic asymptotic description in the {x, t) plane, analytic connection 
formulae have to be established between the various asymptotic formulae. For example, despite 



the asymptotic formulas (2.1) and (2.24) having an error of order e, the numerical results show 



that there is still a region near the leading edge at the boundary of the Whitham zone where the 



error is big ger. T his means that a connection formula between the elliptic expansion (2.1) and the 

I 1 2 ^ 

expansion (2.24) (where the terms of order es have been dropped), is needed. 

In p.3j Claeys has derived connection formula for the P| solution U{X,T) of (2.10) in different 
regions of the {X, T) plane. Using these relations, one can derive in a non rigorous way the 
corresponding connection formulas for the asymptotic solution of the KdV equation near the point 
of gradient catastrophe. Indeed the Pj solution U{X, T) describes a singular transition between 
a region of simple algebraic asymptotics and a region of more complicated oscillatory asymptotics 
involving the Jacobi elliptic ^-function. One may thus expect that U{X, T) itself also exhibits 
different types of asymptotics. Indeed the following result holds |13j : 



X — )■ iboo and T — )• — oo or T — )• +oo in such a way S 



X 

¥1 



remains bounded away from 



4V15 



the interval [— 12-v/3, — - — ], then U{X,T) has an algebraic asymptotics; 
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Figure 20. The difference between the numerical solution to the KdV equation for 
the initial datum uq = — sech'^x at t = tr and the corresponding multiscale solution 
for four values of e, in blue the terms in (2.15) up to order e^/^, in green up to order 



e^/^. Note the scaling of the x and A axes with a factor e^/'' and e^/"^ respectively 
to take care of the expected scaling in x and of the shown error next to the critical 
point. 



if —12^/3 < S < — - — , then U{X,T) has an elliptic asymptotics; 
for S — —12\/3, U{X,T) has a P// asymptotics; 



4vl5 

• for S —7- — - — , U{X,T) has a soliton-like asymptotics. 

Substituting the algebraic asymptotic [13] of the Pj solution into the asymptotic expansion 
(2.15), one obtains in a non rigorous way the connection formula between the Hopf asymptotic 
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(A) 



(B) 



Figure 21. In Fig. (a) the difference between the numerical solution to the KdV 
equation for the initial datum uq = — sech^x and e = 10~^ at t = tc and the 
corresponding P| asymptotic solution in blue, and the difference between KdV and 
Hopf solution in green. In Fig. (b) the edges of the zone where the P| asymptotic 
solution provides a better asymptotic description of KdV than the Hopf solution in 
dependence of e. 



solution (1.4) and (2.15). More precisely in the limit when 



^ _ x-Xc- Qucjt - tc) ^ 



goes to infinity in such a way that s = — o- 

T2 



t-tc 
(A;)3/7el 



Vk ^4 — is outside the interval 



-12a/3. 



4^/15^ 
9 ' 



{t - tc) 



or r — )• — oo then the solution of the KdV equation is approximated by 




where z{s) solves the equation s = 6z — z"^ . The second term on the right hand side of the above 
expression coincides with the solution of the Hopf equation with initial data fiiu) = —ku^. 



8.1. P| asymptotics and elliptic asymptotics. The asymptotic expansions (2.1) and (2.15) 
have a connection region that follows from substituting into (2.15) the elliptic asymptotics for the 
solution obtained in [T3]. In the limit when 



X -Xc- 6uc{t - tc) 
A — ^ , 1 



(A:)i/7 



67 



t-tc 

(A;)3/7el 
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(A) (B) 

Figure 22. In the upper part of Fig. (a) one can see the difference between the 
numerical solution to the KdV equation for the initial datum no = — sech^x and 
e = 10~^ at t = tc and the corresponding Hopf solution. The lower part shows 
the same difference, which is replaced close to the critical point by the difference 
between KdV solution and the P| asymptotic solution (shown in red where the error 
is smaller than the one shown above) . Fig. (b) shows the same situation as the lower 
figure in (a) for two values of e. The A-axis is rescaled by a factor e^/'^, the x-axis 
by a factor e®/"^. 



, , X nrX — Xc — 6uJt — tc) 4vl5 
goes to mnnity m such a way that — 12v3 < s = — 3- = v k 3 < , then the 

T-2 (t-tc)5 9 

solution of the KdV equation is approximated by 



4 



^.1) uix, t, 6) = + (61 + h + bs + 2a + |^(log^)"(j^(s, t); T)) + O ( 



where a, s and T are defined in ( |2.3| ) (with the substitution /3i — )• 6j) and the argument of the 
Jacobi elliptic function ■(? is given by 



^.2) n = _ 2(61 + 62 + 63) - q). 



Here the quantities bi = bi{s), with bi{s) > 62(5) > 63(5), describe the Gurevich-Pitaevski |43J 
self-similar solution of the Whitham equations (2.7) with cubic initial data, namely the function 
q = q{bi,b2,bs) in (2.4) is such that q{b,b,b) = h'^ and the hodograph transform (2.8) takes the 
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Figure 23. The figure shows from top to bottom the numerical solution to the 
KdV equation for the initial datum uq = — sech^a; and e = 10~^ at t = 0.23 > tc, 
the asymptotic solution in terms of Hopf and one-phase KdV solution, the multiscale 



solution (2.24), and the multiscale solution (2.15). 



equivalent form 



6 = -E^'? = ^ [(^1 + ^2 + 63)2 + 2(6? + + 
i=i * 



s = ^(26i -bi + b2 + h)^q + q=^ [ih + 62 + hf - 4(6? + + 61)] 

i=l ^ 



r V(e-6i)(e-62)(^-63)(C + J(6i + 62 + = 0. 



The asymptotic expansion (8.1) should give the connection formula near the point of gradient 
catastrophe (xc,tc) between the one-phase KdV asymptotics (|2.l[) and Pj asymptotics (2.15). It 



is straightforward to check that formula (2.1) reduces to ( |8.ip by expanding the initial data for 
the Whitham equations at the point of gradient catastrophe and keeping the first order correction 
(cubic term). Namely, if i = 1, 2, 3 is the solution of the Whitham equation with the initial data 
( |1.1[ ) and bi is the self-similar solution defined in (8.3), then 



(3i{x, t) = uc + \r-^hi{s) + 0{t - tc). 
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Figure 24. Fig. (a) shows the numerical solution to the KdV equation for the initial 
datum no = — sech'^j; and e = 10~^ at t = 0.23 > tc in blue, the asymptotic solution 



in terms of Hopf and one-phase KdV solution in red, the multiscale solution (2.24) 



in cyan, and the multiscale solution (2.15) in green. Fig. (b) shows the difference 
between the asymptotic solution in terms of Hopf and one-phase KdV solution and 
the numerical solution to the KdV equation in green, between the P// asymptotic 
solution (2.24) and KdV in red, and the Pj asymptotic solution (2.15) and KdV in 
blue. 



With this formula 
limit is of order t — 



1.1) can be recovered from (2.1) by the above substitution. The error in this 
The connection formula 



. 1 ) has already appeared in jS^ . 



8.2. Connection between P^ asymptotic and Pjj asymptotic or soliton asymptotic. 

When S = — 12\/3 the Pj solution has an asymptotic expansion that is provided by oscillations 
whose envelope is given by the Hasting McLeod solution of the Painleve H equation |13] . Plug- 



ging this expansion into (2.15), one obtai ns in a non rigorous way the connection formula between 

I I I I 2 

(2.15) and the P// asymptotic solution (2.24), where the terms of order es have been dropped. 



Introducing the variable 



X + 12^/3T5 _ x-Xc-6uc{t-tc) + 12y/3k~^t-tc)^ _ 2s3T2 



CoCiTs 



esco{t - tc)3 



1 ) 

56 



in the P| asymptotic solution (2.15) and letting T — )■ +oo in such a way that XT- 2 = l2^/3 one 
obtains from [T3] 



^.4) u{x,t,e) = UC + 2V3 



t-tr_ 
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Figure 25. The figure shows from top to bottom the numerical solution to the 
KdV equation for the initial datum uq = —sech^x and e = 10~^ at t = 0.23 > tc, 
the asymptotic solut ion in terms of Hopf and one-phase Kd V so lution, the soliton 
asymptotic solution (2.32), and the asymptotic solution (2.15). 



where q{(,) is the Hasting-McLeod solution of Painleve II, and the phase u is given by 

't-tr_ 



7 



1 ^ - 
k^eT 



This expansion coincides with (2.24) when the initial data at time t = tc is approximated by the 
cubic initial data = —k{u — Uc)^ + 0{u — Uc)"^ with k defined in (2.14). Indeed in this case 

the solution of system (2.19)-(2.21 ) takes the form 



x~{t) -Xc- 6uc{t - tc) = 12\/3^^-^^ + 0{{t - tc 

yk 

+ 0{t-tc), v{t)-Uc = - 



u{t) — Uc = 2y/3\l — j—^ 



\/3 It - tc 



k 



+ Oit-tc). 



Then plugging the above expressions of x (t), u{t) and v{t) into (2.24) one arrives at (8.4) (with 
a different error term though). 

4 



The P| solution has a connection region with the soliton asymptotics (2.32 ) when S = g'V^lS |13j . 

Substituting the corresponding connection formula in {13j into (2.15) one obtains in a non rigorous 
way the connection formula for the asymptotic expansions (2.15) and (2.32). Indeed introducing 
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Figure 26. Fig. (a) shows the numerical solution to the KdV equation for the initial 
datum no = — sech^j; and e = 10~^ at t = 0.23 > tc in blue, the asymptotic solution 



in terms of Hopf and one-phase KdV solution in green, the multiscale solution (2.32 ) 



in red, and the multiscale solution (2.15) in cyan. Fig. (b) shows the difference 
between the asymptotic solution in terms of Hopf and one-phase KdV solution and 
the numerical solution to the KdV in blue, between the soliton asymptotic solution 
(2.32) and KdV in green, and the Pj asymptotic solution (2.15) and KdV in red. 



the variable 



o X - -Vi5tI 

" T^ilogT 



4 3 

6uc{t-tc) - -\/l5fe(i - tc)2 
y 



, N _ 1 , f t — t, 

e{t - tc) 3 log 



£7 A;7 



Co 



(15)4, 



_3 4 

where X and T are defined as in (2.13), and letting T — )• -|-oo in such a way that XT 2 = gVl5 
then 



u{x, t, e) = Uc — 2 



f) ft- tr\ 2 



(8.5) 



3 V A; 



+ 2d 




^sech^Xj + 
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where 



^.6) 



~{]^-i + j)\nT-HV2^h, 



{j + -)log 



16c| 
(15)^ 



2i 



IT 'i^/J\ 



This connection formula coincides with (2.32) when the initial data at i = tc is approximated by 
the cubic initial data fhiu) = —k( u — Uc )"^ + 0{u — Uc)^ where k is defined in (2.14). In this case the 

solution of the system of equation (2.28)- (2.30) takes the form x~^{t) — Xc — Guc{t — tc) = — 



tc)^+0{{t-tc)^), u{t)- 



t-tc 
k 



+ 0{t — tc) and v{t) —Uc 



1 



t-U 



k 



9Vk 

+ 0{t-tc 



Plugging the above expressions of x^{t), u{t) and v(t) into (2.32 ) one obtains the connection formula 

The rigorous derivation of the above connection formulas and their numerical implementation 
will be investigated in a subsequent publication. 
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